R tiene una inmensa capacidad de graficar y visualizar datos de todo tipo, incluídos datos genéticos.

Las gráficas pueden hacerse desde la base de R (base) o con paquetes especializados en graficar, como lattice, o más recientemente ggplot2 y ggbio. También paquetes especializados en un tipo de datos que incluyen funciones para graficar, como ape para árboles filogenéticos.

En esta sección veremos una introducción a graficar en R usando graphics que es el sistema que viene con base y luego nos enfocaremos en gráficas más complejas y las principales usadas en análisis genéticos. En R se puede hacer mucho más que lo que veremos aquí, recomiendo profundizar.

Una de las mejores formas de aprender a hacer gráficas en R es buscar en internet/libro una gráfica parecida a la que queremos hacer y ver el código. Algunas recomendaciones:

Base Graphics

Estas son las principales funciones para graficar utilizando la base de R. Puedes buscar ayuda de cada una con su nombre, y además en explorar argumentos extras con ?par

Gráficas x,y:**

Dando x, y:

largo<-c(10,20,11,15,16,20)
ancho<-c(5,10,7,8,8,11)
plot(x=largo, y=ancho)

Dando un objeto que tiene dos columnas, se toman automático como x,y:

# ver el contenido de `cars`(una df ejemplo que viene con R)
head(cars)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10
plot(cars)

Si queremos especificar qué columnas serán x, y del objeto:

# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist)

Cambiar título de ejes e íconos:

# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist, xlab="Velocidad", ylab="Distancia", cex=0.5, pch=19)

Ejercicio: mira la ayuda de par y explica qué hacen los argumentos cex y pch.

Ejercicio: Repite la figura anterior pero cambiando los puntos por triángulos azules. Necesitarás esto.

Histogramas

Ejemplo con los datos islands (viene con R)

hist(islands)

Ejercicio En “Practicas/Uni7/data/reads.txt” encontrarás un archivo con la cantidad de lecturas de las muestras de tres librerías que fueron secuenciadas en Illumina. Grafica un histograma de las lecturas de cada muestra.

Barplot

Ejemplo:

DNAcon<-data.frame(muestra=c("A", "B", "C"), concentracionADN=c(5,10,9))
barplot(DNAcon$concentracionADN, names.arg=DNAcon$muestra)

Ejercicio Repite la gráfica de DNAcon pero agregando títulos a los ejes x,y

Ahora cargemos un archivo real de datos:

reads<-read.delim("../Practicas/Uni7/data/reads.txt")

Hagamos una gráfica de barras y colorear acorde a info contenida en otra columna:

head(reads)
##   Library    sample  nreads
## 1    Lib1  pobA21_r 1381230
## 2    Lib1    pobB10 1726622
## 3    Lib1    pobB05  819766
## 4    Lib1    pobC05  442508
## 5    Lib1  pobD21_r 1398874
## 6    Lib1 Outa112_r 2024684
barplot(reads$nreads, col=reads$Library)

Definir colores

Los colores que R ocupa para colorear algo están definidos en palette y pueden cambiarse

# Ver colores
palette()
## [1] "black"   "red"     "green3"  "blue"    "cyan"    "magenta" "yellow" 
## [8] "gray"
# Cambiar colores 
palette(c("green", "blue", "red"))

# volver a graficar
barplot(reads$nreads, col=reads$Library)

Además de manualmente, los colores se pueden definir via paletas predeterminadas:

# Cambiar palette a 6 colores del arcoiris
palette(rainbow(6))

# volver a graficar
barplot(reads$nreads, col=reads$Library)

Checa otras palettes parecidas a rainbow en este link, y no te pierdas cómo nombrar muchos otros colores y utilizar otras paletas con más colores en la R Color Reference Sheet

Agregar una leyenda

# Graficar
barplot(reads$nreads, col=reads$Library)
# Agregar leyenda
legend(x="topleft", legend=levels(reads$Library), fill=palette()[1:3])

Nota que legend es una función por si misma (i.e. NO un argumento de plot) que requiere que antes de correrlo se haya corrido plot. Es decir una vez que creamos una gráfica podemos agregar sobre de esta una leyenda. Lo mismo puede hacerse con la función title.

Boxplot

Ejemplo:

boxplot(reads$nreads ~ reads$Library,
        border = c("red", "blue", "darkgreen"))

Ejercicio agrega un título a la gráfica de boxplot de los reads.

Graficar árboles filogenéticos

La graficación de árboles filogenéticos se hace con el paquete ape, con el paquete phytools para funcionalidad extendida y con el paquete ggtree. Empezaremos por ape.

Bibliografía recomendada:

Los árboles filogenéticos pueden construirse en R, simularse en R o leerse a R.

Veamos un ejemplo con un árbol simulado:

# Cargar librería
library(ape)
## Warning: package 'ape' was built under R version 3.2.3
# Simular árbol
set.seed(1) # este comando es opcional, sirve para que todas "simulemos los mismos números" y podamos repetir de forma idéntica la simulación cada vez
tree <- rtree(n = 10, rooted=FALSE)

# ¿Qué tipo de objeto es?
class(tree)
## [1] "phylo"
# ¿Qué contiene?
tree
## 
## Phylogenetic tree with 10 tips and 8 internal nodes.
## 
## Tip labels:
##  t10, t6, t9, t1, t2, t7, ...
## 
## Unrooted; includes branch lengths.
str(tree)
## List of 4
##  $ edge       : int [1:17, 1:2] 11 12 13 13 12 11 14 15 15 14 ...
##  $ tip.label  : chr [1:10] "t10" "t6" "t9" "t1" ...
##  $ edge.length: num [1:17] 0.718 0.992 0.38 0.777 0.935 ...
##  $ Nnode      : int 8
##  - attr(*, "class")= chr "phylo"
##  - attr(*, "order")= chr "cladewise"
# Graficar el árbol
plot.phylo(tree, edge.width=2)

La función plot.phylo puede abreviarse como plot. R sabe que debe usar la función plot.phylo y no plot básico (como lo usamos arriba) porque el objeto que le damos es un árbol.

Podemos modificar este árbol de manera similar a cómo lo hicimos en las gráficas anteriores.

Cambiar el tipo de árbol

Ejercicio Revisa la ayuda de plot.phylo y utiliza un argumento de dicha función para graficar el árbol simulado pero que se vea como un abanico y luego como un cladograma (así):

Enraizar el árbol

# plot árbol sin enraizar
plot.phylo(tree, edge.width=2)

# especificar output para enraizar:
tree<-root(tree, outgroup="t2")

# plot árbol enraizado
plot.phylo(tree, edge.width=2)

Chulear el árbol

Ejercicio: Sigue los ejemplos de Phylogenetic trees in R, del blog Sensory Evolution para realizar los siguientes cambios al árbol anterior:

  • Cambia el nombre de las puntas de “t1”, “t2” etc a “especie 1”, “especie 2”, etc y asegúrate de graficar estos nombres en las puntas de tu árbol
  • Agrega un círculo en los nombres de las puntas. El clado de las puntas t9, t6 y t 10 debe tener círculos rosas, la punta t2 gris y el resto verdes.
  • Incremente al grosor de la línea
  • Cambia el color de la línea a verde oscuro

Leer un árbol en R y graficarlo

Podemos leer árboles en formato newick o nexus a R con la función read.tree de ape:

# cargar archivo
maiz.tree<-read.nexus("../Practicas/Uni7/data/tree")

# checar contenido
maiz.tree
## 
## Phylogenetic tree with 165 tips and 163 internal nodes.
## 
## Tip labels:
##  maiz_3, maiz_68, maiz_91, maiz_39, maiz_12, maiz_41, ...
## 
## Unrooted; includes branch lengths.
# graficar
plot(maiz.tree, type="unrooted", edge.width=0.1, cex=0.5)

Vamos a poner colores de acuerdo a la Categoría de Altitud en vez de nombres de muestras.

### Graficar por Categorías Altitud

# leer info extra de las muestras
fullmat<-read.delim("../Practicas/Uni7/meta/maizteocintle_SNP50k_meta_extended.txt")

# ¿Cuántos colores necesito?
numcolsneeded<-length(levels(fullmat$Categ.Altitud))
palette(rainbow(numcolsneeded)) 

# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)

# Agregar tip labels que correspondan a las categorías de altitud
tiplabels(pch=20, col=fullmat$Categ.Altitud)
# legend
legend(x= "bottomleft", legend=levels(fullmat$Categ.Altitud), pch=19, col=1:numcolsneeded, cex=1, , bty="n")

Ejercicio: Colorea el árbol anterior por Raza (sin incluir una leyenda porque son demasiadas)

Observa que las muestras 162-165 corresponden a teocintles (Zea m. mexicana y Z. m. parviglumis)

fullmat[162:165, 1:4]
##     OrderColecta NSiembra                     Origen               Raza
## 162           NA       96 Jardin Etnobotanico Oaxaca Zea m. parviglumis
## 163           NA      203              Luis Eguiarte    Zea m. mexicana
## 164           NA      911          Luis Eguiarte 1_1 Zea m. parviglumis
## 165           NA     9107          Luis Eguarte 10_7    Zea m. mexicana

Si quisiéramos colorear todas las tips de maíz, pero no las teocintle podemos especificar esto así:

# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)
# Agregar tip labels sólo a los maíces
tiplabels(tip=c(1:161), pch=20, col="black")

Ejercicio Grafica el árbol de maíces de manera que los teocintles sean cuadrados negros y lo s maíces círculos verdes, así:

Mapas en R

En R pueden visualizarse mapas de muchas maneras y de hecho hasta hacer análisis complejos con datos raster, como simulaciones de cambio climático y modelos de distribución potencial de las especies. Aquí sólo cubriremos brevemente cómo graficar un shapefile y agregar puntos.

Carguemos uno de los principales paquetes para manipular mapas:

library(maptools)
## Warning: package 'maptools' was built under R version 3.2.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.2.5
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()

Leer un shapefile

La función readShapePoly de dicho paquete nos permite leer un shapefile de polígonos (para puntos hay que usar otra función parecida ¿cómo crees que se llame?).

Por ejemplo vamos leer a R y graficar el shaphile de las regiones biogeográficas de México:

# leer shapefile
biogeo<-readShapePoly("../Practicas/Uni7/data/rbiog4mgw/rbiog4mgw.shp")

# plot
plot(biogeo)

# colorear por bioregioón
palette("default")
plot(biogeo, border="grey", col=biogeo$PROVINCIA)

## Cambiar colores default a algo más bonito
#¿Cuántos colores necesito?
levels(biogeo$PROVINCIA)
##  [1] "Altiplano Norte (Chihuahuense)"     
##  [2] "Altiplano Sur (Zacatecano-Potosino)"
##  [3] "Baja California"                    
##  [4] "California"                         
##  [5] "Costa del Pacifico"                 
##  [6] "Del Cabo"                           
##  [7] "Depresion del Balsas"               
##  [8] "Eje Volcanico"                      
##  [9] "Golfo de Mexico"                    
## [10] "Los Altos de Chiapas"               
## [11] "Oaxaca"                             
## [12] "Peten"                              
## [13] "Sierra Madre del Sur"               
## [14] "Sierra Madre Occidental"            
## [15] "Sierra Madre Oriental"              
## [16] "Soconusco"                          
## [17] "Sonorense"                          
## [18] "Tamaulipeca"                        
## [19] "Yucatan"
# Generar paleta con colores de RColorBrewer
# ver opciones de colores
library(RColorBrewer)
display.brewer.all()  

# generar paleta
palette(c(brewer.pal(9, "Set1"), brewer.pal(10, "Set3")))

# plot
plot(biogeo, border="grey", col=biogeo$PROVINCIA)
legend("bottomleft", legend=levels(biogeo$PROVINCIA), bty="n", cex=.4, fill=palette())

Ejercicio: cambia el color de las provincias “Tamaulipecas” y “Costa del Pacífico” a otro color.

Agregar puntos a un mapa

Es muy común tener las coordenadas x,y de nuestros puntos de muestreo en un archivo junto con el resto de la info de nuestras muestras. Por ejemplo en el caso de la info de maíces que hemos estado utilizando:

fullmat[1:5,c("Latitud", "Longitud")]
##    Latitud   Longitud
## 1 28.68578 -107.58300
## 2 19.29083  -97.92944
## 3 28.53703 -108.92467
## 4 16.03225  -92.97972
## 5 20.09111 -101.11889

Agregar esta info a un mapa en la misma proyección y sistema de coordenadas puede hacerse con la función points:

# plot map
plot(biogeo, border="grey", col=biogeo$PROVINCIA, lwd=0.8)
# agregar puntos
points(fullmat$Longitud, fullmat$Latitud, pch=19, col="black", cex=0.4)

Ejercicio: Baja un mapa (nivel nacional) que te interese del GeoPortal de la CONABIO, ploteálo y agrega los puntos del muestro de maíz, utilizando una forma de punto diferente para los teocintles, y que los puntos estén coloreados por CategAltitud.

Utilizar mapas de GoogleMaps

Hay una librería para bajar mapas (estáticos) de GoogleMaps: RgoogleMaps

Ejemplo:

# Cargar librería
library(RgoogleMaps)

# Establecer limites deseados
lat <- c(13,30) # definir limites en y
lon <- c(-112,-93) #definir limites en x
center = c(mean(lat), mean(lon))  #donde queremos centrar
zoom <- 4  #zoom deseado (1 lo más alejado)

# obtener mapa
mygooglemap <- GetMap(center=center, zoom=zoom, maptype= "terrain")

# plot
PlotOnStaticMap(mygooglemap)

Guardar imágenes:

Método 1:

  • “Abrir un device” con png(), jpg(), pdf() según el formato en que queramos guardar.
  • Hacer la gráfica
  • “Cerrar el device” con dev.off
png("../Practicas/Uni7/out/arbol.png")
plot(maiz.tree, type="unrooted", show.tip=FALSE)
tiplabels(pch=20, col="green3", tip=c(1:161), cex=.5)
tiplabels(pch=15, col="black", tip=c(162:165), cex=.5)
dev.off()
## quartz_off_screen 
##                 2

Método 2:

En R studio darle “Export” en el panel de la imagen y seleccionar un nombre de archivo y demás características.

ggplot2

Las gráficas que hemos visto hasta ahora pueden verse un poco feas de inicio y puede tomar un rato y mucho código arreglarlas a algo hermoso. ggplot2 es un paquete que ahorra este trabajo y que ha comenzado a ser ampliamente adoptado.

ggplot2 construye gráficas “definiendo sus componentes paso a paso”.

Para poder usar ggplot2 se requiere que la data.frame esté en formato largo, como vimos cuando revisamos la función gather. Además de esos apuntes puedes revisar esto si te quedan dudas.

Términos importantes:

Ojo: Mucho mejor que ver la ayuda de cada función es ver la Documentación online de ggplot2 y este R Graphics Cookbook

Ejemplos:

Gráficas de dispersión

# Cargar ggplot2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
# Examinar datos pre-cargados
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
# graficar 
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) + geom_point()

Pregunta: ¿Qué hace el símbolo +? Nota que el código anterior tmb puede escribirse así:

myplot<-ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width))

myplot + geom_point()

Los colores y formas se cambian en aes:

ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width, color= Species, shape=Species)) + geom_point()

Ya sea en el aes de la función inicial o dentro de los geoms (Nota que el tamaño no es un aes, sino un argumento de geom_point)

ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) + 
    geom_point(aes(color= Species, shape=Species), size=3)

Si queremos quitar el fondo gris:

ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) + 
    geom_point(aes(color= Species, shape=Species), size=3) + 
    theme_bw()

Aveces queremos graficar en páneles separados la misma info para diferentes tratamientos o especies. Por ejemplo:

ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) + 
  geom_point() +
  facet_grid(Species ~ .)

Ejercicio Pon color por especie a la gráfica anterior:

Ejercicio Repite la gráfica anterior pero para que se vea así:

Ejercicio Repite la figura anterior pero cambiando los labels para que digan “Ancho de sépalo” y “Largo de sépalo”, respectivamente. Debe verse así:

También podemos agregar el resultado de un modelo matemático, como una regresión lineal:

Barplot

En este tipo de gráficas la altura de las barras puede significar dos cosas:

  • la cuenta (frecuencia) de casos de cada valor de x. Si quieres graficar esto utiliza stat="bin" en los argumentos de geom_bar.
  • el valor de la columna en el set de datos. Si quieres graficar esto utiliza stat="identity" en los argumentos de geom_bar.

Ejemplo:

# Cargar archivo
reads<-read.delim("../Practicas/Uni7/data/reads.txt")

# plot
p <- ggplot(data=reads, aes(x=sample, y=nreads, fill=Library)) + 
  geom_bar(stat="identity") 
p 

# Rotar nombres muestras
p <- p +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=6))
p

La gráfica anterior no es igual a la que hicmos con barplot con los mismos datos ya que ggplot2 grafica en el orden de los levels, en este caso:

head(levels(reads$sample))
## [1] "NegC1"     "NegC2"     "NegC3"     "Outa112"   "Outa112_r" "Outa113"

Forma de solucionarlo:

# Cambiar orden de levels:
reads$sample<-factor(reads$sample, levels = reads$sample[order(1:length(reads$sample))])
head(levels(reads$sample))
## [1] "pobA21_r"  "pobB10"    "pobB05"    "pobC05"    "pobD21_r"  "Outa112_r"
# Graficar
# plot
p <- ggplot(data=reads, aes(x=sample, y=nreads, fill=Library)) + geom_bar(stat="identity") 
p +  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=6))

Boxplot

Siguiendo con los mismos datos anteriores:

# plot
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p

# quitar leyenda
p + guides(fill=FALSE)

Árboles

Los árboles filogenéticos también pueden graficarse estilo ggplot, con el paquete de Biocounductor ggtree, que veremos más adelante junto con otros paquetes.

Cambiar colores en ggplot

Al igual que en base, en ggplot es posible cambiar los colores manualmente o cambiando la paleta.

Recomiendo buscar más información y ejemplos en esta excelente guía Cookbook-R colores en ggplot2.

Ejemplos:

Cambiar colores manualmente

# Crear paleta:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Usar paleta en gráfica:
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_manual(values=c("red", "blue", "green"))

Cambiar la paleta

# Crear paleta apta para daltónicos:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Usar paleta en gráfica:
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_manual(values=cbPalette)

Usar una paleta de ColorBrewer

p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_brewer(palette="Set1")

Utilizar gradientes de colores

# Generar datos
set.seed(1)
df <- data.frame(xval=rnorm(50), yval=rnorm(50))

# Plot
ggplot(df, aes(x=xval, y=yval, colour=yval)) + geom_point()

# Cambiar gradiente
ggplot(df, aes(x=xval, y=yval, colour=yval)) + geom_point() + 
    scale_colour_gradientn(colours=rainbow(6))